Practical image reconstruction for magnetic resonance imaging

ABSTRACT

A set of image-space data is reconstructed from a set of k-space data. The set of image-space data is generated by minimizing a cost functional by an iterative non-linear conjugate gradient process. The iterative process may be accelerated by introducing k-space weighting to the cost functional. With proper choice of k-space weighting, a block-Toeplitz matrix is generated which permits use of Fast Fourier Transform techniques. An image is rendered from the set of image-space data.

This application claims the benefit of U.S. Provisional Application No.60/864,189 filed Nov. 3, 2006, which is incorporated herein byreference.

BACKGROUND OF THE INVENTION

The present invention relates generally to medical imaging, and moreparticularly to iterative reconstruction of images.

Magnetic Resonance Imaging (MARI) has become a well-established medicaldiagnostic tool for imaging structures within the body of a patient.Image quality may be characterized by a host of parameters, includingresolution, field of view, contrast, edge definition, and artifacts (forexample, ghosts and streaks). Under a broad range of conditions, imagequality improves with increasing data acquisition time. If the dataacquisition time is increased, however, the patient is subjected to alonger scan time, which increases patient discomfort. In some instances,long scan times may actually degrade image quality because of movementof the region of interest during the scan. Short scan times are alsonecessary for near-real-time measurements, such as used in functionalMRI. There is, thus, a trade-off between image quality and scan time.

Images are displayed on physical media; for example, emulsions on filmor pixels on a monitor. The normal physical world may be referred to asreal space. In one method for producing high-quality images, MR signalsare captured in k-space. In some fields of study, k-space is alsoreferred to as spatial-frequency domain. In general terms, data valuesin real space are then generated by taking the inverse Fourier transformof data values in k-space. In general, MR signals are not measured as acontinuous function of position in k-space. They are sampled at discretek-values. Subject to specific constraints and boundary conditions, imagequality generally improves as the density and range of discrete k-spacesampling points are increased. Recording a large number of samples,however, has disadvantages. One is the extended scan time discussedabove. The other is low temporal resolution.

To reduce data acquisition time, MRI data are often intentionallyunder-sampled. This will, however, often lead to reduced signal-to-noiseratio (SNR) and to image degradation. Various techniques have beendeveloped to enhance the image quality reconstructed from under-sampleddata, but they require extended computational time and high memoryusage. What is needed is a method which reduces computational time andmemory usage for generating high quality real-space images fromunder-sampled k-space data.

BRIEF SUMMARY OF THE INVENTION

An image is rendered from a set of k-space data by first reconstructinga set of image-space data from the set of k-space data. The set ofimage-space data is generated by iteratively solving for the minimum ofa cost functional comprising a weighting matrix. With proper choice ofthe weighting matrix, a block-Toeplitz matrix is generated. The set ofk-space data and an initial image-space data are inputted into anon-linear conjugate gradient process for iteratively reconstructingimage-space data. The k-space weighting matrix and the block-Toeplitzmatrix accelerate the iterative solution process. These and otheradvantages of the invention will be apparent to those of ordinary skillin the art by reference to the following detailed description and theaccompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 shows a high-level schematic of an MRI system;

FIG. 2 shows a high-level flowchart of an image reconstruction process;and,

FIG. 3. shows a high-level schematic of a computer for performing imagereconstruction.

DETAILED DESCRIPTION

FIG. 1 shows a high-level schematic of the functional architecture of arepresentative imaging system. In the example shown in FIG. 1, theimaging system is MRI system 100. Embodiments apply to other modalities.MRI system 100 comprises image measurement system 102, image processingsystem 110, and display 118. Image measurement system 102 may includescanning system 104 and signal processing system 106. Scanning system104 includes the MRI scanner and associated instrumentation andsoftware. The output of scanning system 104 are radio-frequency (RF)signals which are fed into signal processing system 106. Signalprocessing system 106 includes hardware and software which converts theRF signals to k-space data 108, which is discussed further below.

The k-space data 108 is then fed into image processing system 110, whichmay include image reconstruction system 112, image rendering system 114,and video processing system 116. Image reconstruction system 112,further discussed below, transforms k-space data 108 into image-spacedata, which may be a two-dimensional (2-D) or three-dimensional (3-D)dataset. Herein, image space refers to real space. Herein, imagereconstruction refers to the transformation of k-space data toimage-space data. The image-space data is then mapped by image renderingsystem 114 into optical data for further video processing. For example,in a monochrome display, the image-space data may be mapped into theluminance values of pixels. In a color display, image-space data may bemapped into the luminance and false-color values of pixels. Videoprocessing system 116 transforms the output optical data from imagerendering system 114 into signals which drive pixels on display 118,which displays the image viewed by a user.

As discussed above, MR signals captured and processed by imagemeasurement system 102 are output as k-space data 108. In general terms,data values in image space are then generated by taking the inverseFourier transform of k-space data 108. In general, MR signals are notmeasured as a continuous function of position in k-space. They aresampled at discrete k-values. Subject to specific constraints andboundary conditions, image quality generally improves as the density andrange of discrete k-space sample points are increased. For othermodalities, k-space is also referred to as spatial-frequency domain.

Coordinates of points in k-space may be specified with respect to a 3-DCartesian grid, analogous to the common x-y-z Cartesian grid in realspace. For MRI analysis, the axes are k_(x)=frequency (also calledreadout), k_(y)=phase encoding, and k_(z)=slice selection. Alternativecoordinate systems, such as polar coordinates, may also be used. Thelocus of points in k-space along which data is sampled is referred to asthe scan trajectory. The simplest trajectories are linear. For example,in a 2-D scan, a value of k_(y) is held fixed while data is sampled atdifferent values of k_(x). The linear scans along k_(x) are thenrepeated for other fixed values of k_(y). A disadvantage of linearCartesian scans is that the scan pattern is relatively slow in the k_(y)direction.

Acquisition speed may be increased by various techniques. For example,radial or spiral trajectories may be used instead of linear ones. Asdiscussed above, k-space values may be intentionally under-sampled.Images reconstructed from under-sampled data, however, are oftendegraded by severe artifacts and low SNR. Image quality may be enhancedvia iterative MR reconstruction schemes. (See, for example, T-C Chang,et al., MR Image Reconstruction from Sparse Radial Samples Using BregmanIteration, Proc. 14th. Annual Meeting of ISMRM, Seattle, Wash., 2006.)In these schemes, the L₁ norm of the sparse representations of theimage-space data is minimized, subject to the constraint of k-space datafidelity. A sparse representation of an image is a set of transformeddomain coefficients that has much fewer non-zero values than the imageof interest. Using the set of coefficients, the image of interest may beconstructed, rendered, and displayed.

The problem can be expressed as:

$\begin{matrix}{\begin{matrix}{\min\limits_{f}{{\varphi(f)}}_{1}} & {{{{subject}\mspace{14mu}{to}\mspace{14mu}{Af}} = y},}\end{matrix}\mspace{14mu}} & \left( {E\; 1} \right)\end{matrix}$where ƒ is the image-space data arranged in a vector form, φ(□)transforms the image-space data into a sparse representation, y is themeasured k-space data arranged in a vector form, and matrix A is aFourier operator that maps the image-space data into k-space data. Here,∥•∥_(p) refers to the L_(p) norm. The image is rendered from ƒ-valueswhich solve the minimization E1.

The optimization problem of E1 is equivalent to the unconstrainedoptimization:

$\begin{matrix}{{\min\limits_{f}\left\{ {\psi(f)} \right\}},{{{with}\mspace{14mu}{\psi(f)}} = {{\lambda{{\varphi(f)}}_{1}} + {{{{Af} - y}}_{2}^{2}.}}}} & \left( {E\; 2} \right)\end{matrix}$Here, ψ(ƒ) is a cost functional, λ≧0 is a Lagrange multiplier, and∥Aƒ−y∥₂ ² is a data fidelity term. In an embodiment, an iterativenon-linear conjugate gradient (CG) method is used to calculate ƒ-valueswhich solve the minimization E2. The solution at iteration k isƒ_(k+1)=ƒ_(k)+α_(k) d _(k),  (E3)whered _(k+1)=−∇ψ(ƒ_(k+1))+β_(k+1) d _(k)  (E4)d ₀=−∇ψ(ƒ₀)  (E5)β_(k+1) is the update direction:

$\begin{matrix}{{\beta_{k + 1} = \frac{{{\nabla\;{\psi\left( f_{k + 1} \right)}}}^{2}}{{{\nabla\;{\psi\left( f_{k} \right)}}}^{2}}},} & \left( {E\; 6} \right)\end{matrix}$α_(k) is the step size:

$\begin{matrix}{\alpha_{k} = {\underset{\alpha}{\arg\;\min}\;{{\psi\left( {f_{k} + {\alpha\; d_{k}}} \right)}.}}} & \left( {E\; 7} \right)\end{matrix}$In E7, a search for α_(k) may be performed with a backtracking linesearch method. An example of a backtracking line search method is givenby the pseudo-code:

given a descent direction d for the energy functional ψ, γε(0, 0.5),δε(0, 1);

initializing α:=1;

and while ψ(ƒ+αd)>ψ(ƒ)+γα∇ψ(ƒ)^(T)d, setting α:=δα.

Using this formulation, the number of iterations necessary forconvergence is high. The cost of each CG iteration is high as well, dueto the non-Cartesian nature of the trajectory. The gradient of the datafidelity term in E2 is A^(†)(Aƒ−y). Fast Fourier transform (FFT) is notdirectly applicable when applying the Fourier operator A to theimage-space data and its adjoint A^(†)to the k-space data.Conventionally, gridding and inverse gridding methods are used. (See,for example, J. I. Jackson, et al., Selection of a Convolution Functionfor Fourier Inversion Using Gridding, IEEE Trans. Med. Imag., 10(3),473-478, 1991.) The gridding accuracy trades off computational speedsubject to the choice of gridding kernel. Furthermore, the use of alook-up table in the gridding methods trades off memory usage forcomputational speed. As a result, using gridding and inverse griddingfor the data fidelity term is expensive for 2-D applications, and may beprohibitively high for 3-D applications.

Image reconstruction time may be greatly reduced by embodiments of theinvention. In an embodiment, a k-space weighting, equivalentlyinterpreted as a density compensation function (DCF), is applied to theL₂ norm data fidelity term in E2. The minimization problem of E2 is thenmodified to

$\begin{matrix}{{\min\limits_{f}\left\{ {{\lambda{{\varphi(f)}}_{1}} + {{W\left( {{Af} - y} \right)}}_{2}^{2}} \right\}},} & \left( {E\; 8} \right)\end{matrix}$where W is a diagonal weighting matrix. The modified cost functional isλ∥φ(ƒ)∥₁+∥W(Aƒ−y)₂ ². The matrix D=W⁺W can be proportional to thedensity compensation function, which is well-known in gridding method.The weighting makes the residual nearly perpendicular to the descentdirection at each iteration, and thus enables near optimal convergence.In an embodiment, the product of the Fourier operator and its adjoint isa block-Toeplitz matrix, and its operation on vectors may be evaluatedby FFT. Gridding and inverse gridding time are reduced. The storage ofthe gridding look-up table is also eliminated.

In an embodiment, total variation (TV) may be used as ∥φ(□)∥₁ due to itsconvexity, simplicity, and effectiveness for processing medical images.The models represented by E2 and E8 then become total variationregularization models, which are advantageous in a wide variety ofrestoration problems, such as denoising, deblurring, blinddeconvolution, and in-painting. TV norm is effective in reducingoscillations without penalizing discontinuities. Thus, it may mitigatestreaking artifacts while preserving edge definition. By selectingdifferent values for λ, a user may select the level of detail desired inthe reconstructed image. The Lagrange multipliers λ is usually chosenrelatively small to emphasize the data fidelity term and avoid artifactsintroduced by over-penalizing the TV norm. Examples of artifacts includestaircasing (blocking) and degradation of contrast, geometry, andtexture. The convergence of the algorithm is, therefore, largely decidedby the data fidelity term.

In an embodiment, as shown above in E8, a diagonal weighting matrix W isincorporated in the data fidelity term. A value for W that achieves nearoptimal convergence may be calculated by the following procedure.Referring to E3 and given the ƒ-values after the (k−1)^(th) iterationƒ_(k−1), the k^(th) image, without considering the TV gradient (thefirst term in E8) and any scalar factor, is approximately

$\begin{matrix}\begin{matrix}{f_{k} \approx {f_{k - 1} + \left( {{{- A^{\dagger}}{DAf}_{k - 1}} + {A^{\dagger}{Dy}}} \right)}} \\{= {{A^{\dagger}{Dy}} + {\left( {I - {A^{\dagger}{DA}}} \right)f_{k - 1}}}}\end{matrix} & \left( {E\; 9} \right)\end{matrix}$

-   -   where I refers to the identity matrix.        Herein, the k^(th) image refers to set of ƒ-values=ƒ_(k).        To achieve optimal convergence, the following residual is set        equal to 0:

$\begin{matrix}\begin{matrix}{{{Af}_{k} - y} = {{{AA}^{\dagger}{Dy}} + {\left( {A - {{AA}^{\dagger}{DA}}} \right)f_{k - 1}} - y}} \\{= {{\left( {{{AA}^{\dagger}D} - I} \right)y} + {\left( {I - {{AA}^{\dagger}D}} \right){Af}_{k - 1}}}} \\{{= {{\left( {{{AA}^{\dagger}D} - I} \right)\left( {{Af}_{k - 1} - y} \right)} = 0}},}\end{matrix} & \left( {E\; 10} \right)\end{matrix}$For Aƒ_(k−1)≠y, the solution of E10 is

$\begin{matrix}{D = \left( {AA}^{\dagger} \right)^{- 1}} & \left( {E\; 11} \right) \\{W = {\left( {AA}^{\dagger} \right)^{\frac{1}{2}}.}} & \left( {E\; 12} \right)\end{matrix}$

The physical meaning of this solution may be clarified by comparing thematrix form of gridding method with a least-squares solution of the MRsignal model (see R. V. de Walle, et al., Reconstruction of MR Imagesfrom Data Acquired on a General Nonregular Grid by PseudoinverseCalculation, IEEE Trans. Med. Imag., 12(19): 1160-1167, 2000).

Based on gridding,I _(g) =A ^(†) D′y, D′=diag(d _(i))  (E13)where I_(g) is the image from gridding reconstruction.Based on least squares,I _(LS) =A ^(†)(AA ^(†))⁻¹ y  (E14)where I_(LS) is the image from least-squares reconstruction.If I_(g) is set equal to I_(LS), the density compensation matrix in E13is a diagonal approximation of (AA^(†))⁻¹ in E14. The same approximationmay also be used to find the weighting matrix D. The weighting matrix Wis then found by taking the square root of D's diagonal elements. Thecalculation of matrix D and matrix W are further discussed below.

Another interpretation of the solution (E11 and E12) is that applyingthe weighting W orthogonalizes the Fourier operator A. For Cartesiansampling, Fourier operator A is orthogonal, meaning that the output oftaking the FFT and then inverse FFT will yield an identical result asthe input. With non-Cartesian sampling (e.g. radial sampling) A is nolonger orthogonal. With E12, however, the weighted Fourier operator B=WAsatisfies the equationB ⁺ B=(WA)⁺ WA=A ^(†) W ⁺ A=A ^(†)(AA ^(†)) ⁻¹ A=I.  (E15)This suggests that a near optimal convergence is achieved when theprojection operator A in the general model, as represented in E1, isorthogonal. Compared to the model represented in E2, the modelrepresented in E8 has a relatively greater emphasis on the highfrequency measurements. This suggests that the new optimization,represented in E8, will favor high resolution fine details rather thanlow resolution contrast in the image. It has the potential, however, foramplifying noise or other artifacts. The Lagrange parameter λ may bechosen relatively large to reduce these amplified oscillations.

While calculating the gradient and the cost functional represented inE8, the Fourier operator A is frequently applied to the image-spacedata, followed by its adjoint A^(†). In the following example, theinstance of radial sampling in k-space is discussed. If a 2-D N×N imageis arranged as a vector of N² length, and if there are N_(r) radiallines, with N₊samples acquired for each radial line, the matrix A thenhas a size of (N_(r)N₊)×N². Direct evaluation of Aƒ in 2-D applicationsis computationally intense. Three-dimensional applications are even moreintensive. One approach to accelerating this evaluation is throughgridding and inverse gridding which approximates A and A^(†)withA=TFP, A^(†)P⁺F⁺T⁺,  (E16)where T is a sparse interpolation matrix with non-zero entries definedby the kernel function

${t_{j,i} = {\Phi\left( {x_{j} - \frac{l}{\mu\; N}} \right)}},$F is a regular Cartesian Fourier matrix which may be efficientlycalculated using FFT, and P is deapodization diagonal matrix assembledby the Fourier transform of Φ. The gridding accuracy is controlled bythe choice of kernel function Φ, and its parameters L (kernel size), andμ (over-sampling rate). As the values of L and μ increase, the accuracyof the approximation increases. The computational time and memory usagehowever, also increase. The matrix T is usually pre-constructed andstored in memory to avoid the kernel calculation on the fly. Although Tis sparse, it still occupies a large memory, especially in 3-D instancesand instances with large values of L.

In an advantageous embodiment, the block-Toeplitz nature of the productA^(†)DA is used to reduce computational time and memory usage. TheToeplitz characteristic of A^(†)A has been previously used for solving ageneral non-uniform inverse FFT (NU-IFFT) problem and otherapplications. (See, for example, Q. H. Liu and X. Y. Tang, IterativeAlgorithm for Nonuniform Inverse Fast Fourier Transform (NU-IFFT),Electronic Letters, 34: 1913-1914, 1998, and J. A. Fessler, et al.,Toeplitz-Based Iterative Reconstruction for MRI with Correction forMagnetic Field Inhomogeneity, IEEE Trans. Sig. Proc., 53: 3393-3402,2005.) The matrix A^(†)DA may be written as

$\begin{matrix}{{A^{\dagger}{DA}} = {\sum\limits_{p,{q \in \Omega}}{\sum\limits_{m = I}^{M}\;{d_{m}{\mathbb{e}}^{{\mathbb{i}2\pi}\;{{k_{m}{({p \cdot q})}}/N}}}}}} & \left( {E\; 17} \right)\end{matrix}$where d_(m) is the m^(th) element on the diagonal of weighting matrix D,k_(m) is the k-space coordinates of the m^(th) sample on the trajectory,and p, q are indexes of image-space Ω. E17 shows that the elements ofA^(†)DA are only dependent on p−q. Therefore, only a small portion ofelements are independent. A sub-matrix assembled by these independentelements is denoted by Γ. The sub-matrix Γ has a matrix size ofM_(Γ)=(2N−1)^(d). Furthermore, the conjugate symmetry givesΓ_(n)=Γ*_(−n). Consequently, only (M_(Γ)+1)/2 elements, out of N^(2d),need to be calculated in a d-dimension case. They can be eitherevaluated exactly using E17 or approximately using gridding method.

Applying the operator A^(†)DA on an image vector ƒ becomes ad-dimensional convolution between the matrix Γ and the d-dimensionalimage ƒ^(d). According to the convolution theorem, it can be computedfast with d-dimensional FFT. Therefore, we have

$\begin{matrix}{{A^{\dagger}{DAf}} = {{\sum\limits_{p \in \Omega}{\sum\limits_{q}{\Gamma_{p - q}f_{q}^{d}}}} = {{{IFFT}\left( {{{FFT}(\Gamma)} \cdot {{FFT}\left( f^{d} \right)}} \right)}.}}} & \left( {E\; 18} \right)\end{matrix}$

Above, IFFT is the inverse FFT operator. Note that A^(†)DA is onlydependent on the sampling trajectory, not the measurements. Thus, FFT(Γ) can be evaluated only once beforehand and then used for all theiterations. Each calculation of E16, therefore, essentially consists ofone FFT and one IFFT. Compared to gridding and inverse griddingapproach, the interpolation and deapodization steps are eliminated, andthe storage of interpolation matrix T is avoided.

FIG. 2 shows a high-level flowchart for an embodiment of a process forreconstructing an image from a set of k-space data. The processcomprises two major phases, a one-time initialization phase 202 and anon-linear conjugate gradient (CG) iteration phase 216. In step 204, agridding table is built in random access memory (RAM). In step 206, aweighting matrix D is calculated. In step 208, a gridding image iscalculated. The gridding image is the density compensated reconstructionimage I_(g)=A^(†)Dƒ. In step 210, the independent elements of theToeplitz matrix A^(†)DA are calculated. After the above values have beencalculated, gridding is no longer needed. Therefore, in step 212, thegridding table is deleted to reduce RAM usage. The values calculatedfrom the initialization phase are provided as input to step 218 in theCG iteration phase 216.

The CG iteration phase 216 comprises a series of iterations of steps218-224. In step 218, using the values calculated in steps 206-210 ofthe initialization phase 202, the gradient ∇ψ=λ∇∥φ(ƒ)∥₁+2A^(†)D(Aƒ−y) iscalculated. In step 220, a backtracking line search process is thenperformed to find a proper step size using the Fletcher-Reeves formula(see, for example, J. E. Shewchuk, An Introduction to the ConjugateGradient Method Without the Agonizing Pain, Edition 1-1/4, School ofComputer Science, Carnegie Mellon University, Aug. 4, 1994, pp. 48 and52). In step 222, the cost functional ψ(ƒ)=λ∥φ(ƒ)∥₁+∥W(Aƒ−y)∥₂ ² iscalculated. In step 224, the image valued is updated. Steps 220-224 arethen iterated until stopping criteria is attained. Stopping criteria,for example, may be a predetermined number of iterations or apredetermined rate of convergence.

Examples of processes for some of the above steps are given herein.

The weighting matrix may be derived using various criteria. For example,an iterative method proposed by Pipe may be used. (See, for example, J.G. Pipe, Sampling Density Compensation in MRI: Rationale and anIterative Numerical Solution, Magn. Reson. Med., 41: 179-186, 1999, andJ. G. Pipe, Reconstructing MR Images from Undersampled Data:Data-Weighting Considerations, Magn. Reson. Med., 43: 867-875, 2000.)Pipe's work indicated that the density compensation function for a 2-Dunder-sampled radial trajectory should be linear within the region thatis sampled above the Nyquist rate, and becomes flat in the under-sampledperiphery region. The DCF is found by recursively computing:

${W_{i + 1} = \frac{W_{i}}{W_{i} \otimes C}},$with an initial

$\begin{matrix}{W_{1} = \frac{S}{S \otimes C}} & \left( {E\; 19} \right)\end{matrix}$where C is a circular convolution kernel (bell shape), and S is thek-space sampling function. The iterations stop when E=S×(W

C) becomes a 1 vector and the range of variation is considerably small,e.g., (0.99, 1.01).

An example of a process for calculating the Toeplitz matrix is describedherein. The Toeplitz matrix is only trajectory dependent, and it iscomputed during the initialization phase. To calculate the (2N−1)^(d)matrix

$\begin{matrix}{{\Gamma_{n} = {\sum\limits_{m = 1}^{M}{d_{m}{\mathbb{e}}^{{\mathbb{i}2\pi}\;{k_{m} \cdot {n/N}}}}}},{n \in \left( {{{- N} + 1},{N - 1}} \right)^{d}},} & \left( {E\; 20} \right)\end{matrix}$it is divided into 2^(d) blocks, half of which can be obtained viaconjugate symmetry from the other half. These blocks are calculated byapplying gridding method to reconstruct a phase-shifted DCF vector. Forexample, the nε(0, N−1)^(d) block is calculated as

$\begin{matrix}{\begin{matrix}{{\Gamma_{n} = {\sum\limits_{m\mspace{11mu} = \mspace{11mu} 1}^{\; M}\;{d_{m}{\mathbb{e}}^{\;{{\mathbb{i}\pi}\; k_{m}}}{\mathbb{e}}^{{\mathbb{i}2\pi}\;{k_{m} \cdot {r/N}}}}}},} & {r \in \left( {{{- N}/2},{{N/2} - 1}} \right)^{d}}\end{matrix}.} & \left( {E\; 21} \right)\end{matrix}$

An example of calculating the gradient of a cost functional is describedherein. In the general 3-D instance, the gradient of the cost functionalis∇ψ=λ∇∥φ(ƒ)∥₁+2A ^(†) D(Aƒ−y).  (E22)It includes the TV term

$\begin{matrix}{{{\nabla{{\varphi(f)}}_{1}} = {\varphi_{{i - 1},j,k}^{x} - \varphi_{i,j,k}^{x} + \varphi_{i,{j - 1},k}^{y} - \varphi_{i,j,k}^{y} + \varphi_{i,j,{k - 1}}^{z} - \varphi_{i,j,k}^{z}}},} & \left( {E\; 23} \right) \\{where} & \; \\\begin{matrix}{\varphi_{i,j,k}^{x} = \frac{\left( {\nabla f_{x}} \right)_{i,j,k}}{\sqrt{\left( {\nabla f_{x}^{2}} \right)_{i,j,k} + \left( {\nabla f_{y}^{2}} \right)_{i,j,k} + \left( {\nabla f_{z}^{2}} \right)_{i,j,k} + \eta}}} \\{\varphi_{i,j,k}^{y} = \frac{\left( {\nabla f_{y}} \right)_{i,j,k}}{\sqrt{\left( {\nabla f_{x}^{2}} \right)_{i,j,k} + \left( {\nabla f_{y}^{2}} \right)_{i,j,k} + \left( {\nabla f_{z}^{2}} \right)_{i,j,k} + \eta}}} \\{\varphi_{i,j,k}^{z} = \frac{\left( {\nabla f_{z}} \right)_{i,j,k}}{\sqrt{\left( {\nabla f_{x}^{2}} \right)_{i,j,k} + \left( {\nabla f_{y}^{2}} \right)_{i,j,k} + \left( {\nabla f_{z}^{2}} \right)_{i,j,k} + \eta}}}\end{matrix} & \left( {E\; 24} \right)\end{matrix}$and the L₂ data fidelity termA ^(†) D(Aƒ−y)=A ^(†) DAƒ−A ^(†) Dy=TP(ƒ)−I _(g),  (E25)where TP(ƒ) stands for the Toeplitz operation, and I_(g) is the densitycompensated gridding reconstruction.

A process for calculating the cost functional is described herein. As ageneral instance, the 3-D cost functional isψ(ƒ)=λ∥φ(ƒ)∥₁ +∥W(Aƒ−y)∥₂ ²  (E26)It again includes the TV term

$\begin{matrix}{{{\varphi(f)}}_{1} = {\sum\limits_{i,j,k}\sqrt{\left( {\nabla f_{x}^{2}} \right)_{i,j,k} + \left( {\nabla f_{y}^{2}} \right)_{i,j,k} + \left( {\nabla f_{z}^{2}} \right)_{i,j,k} + \eta}}} & \left( {E\; 27} \right)\end{matrix}$and the data fidelity term

$\begin{matrix}\begin{matrix}{{{W\;\left( {{Af}\; - \; y} \right)}}_{2}^{2} = {\left( {{Af} - y} \right)^{\dagger}W^{\dagger}\;{W\left( {{Af} - y} \right)}}} \\{= {{f^{\;\dagger}\; A^{\dagger}\;{DAf}}\; - \;{y^{\;\dagger}\;{DAf}}\; - \;{f^{\;\dagger}\; A^{\;\dagger}\;{Dy}}\; + \;{y^{\;\dagger}\;{Dy}}}} \\{= {{f^{\dagger}A^{\dagger}{DAf}} - {2\;{{Re}\left( {f^{\dagger}A^{\dagger}{Dy}} \right)}} + {y^{\dagger}{Dy}}}} \\{= {{f^{\dagger} \cdot {{TP}(f)}} - {2\;{{Re}\left( {f^{\dagger}I_{g}} \right)}} + {y^{\dagger}{Dy}}}}\end{matrix} & \left( {E\; 28} \right)\end{matrix}$where Re(□) indicates the operation for taking the real part. Note thatthe most time-consuming part, TP(ƒ), is shared by both the gradient andcost functional computation, and the left operations are mainly thecalculations of vector norms.

One embodiment of an image processing system 110 (FIG. 1) may beimplemented using a computer. As shown in FIG. 3, computer 302 may beany type of well-known computer comprising a central processing unit(CPU) 306, memory 304, data storage 308, and user input/output interface310. Data storage 308 may comprise a hard drive or non-volatile memory.User input/output interface 310 may comprise a connection to a userinput device 322, such as a keyboard or mouse. As is well known, acomputer operates under control of computer software which defines theoverall operation of the computer and applications. CPU 306 controls theoverall operation of the computer and applications by executing computerprogram instructions which define the overall operation andapplications. The computer program instructions may be stored in datastorage 308 and loaded into memory 304 when execution of the programinstructions is desired. Computer 302 may further comprise a signalinterface 312 and a video display interface 316. Signal interface 312may transform incoming signals, such as from measurement system 102(FIG. 1), to signals capable of being processed by CPU 306. Videodisplay interface 316 may transform signals from CPU 306 to signalswhich may drive video display 318. Computer 302 may further comprise oneor more network interfaces. For example, communications networkinterface 314 may comprise a connection to an Internet Protocol (IP)communications network 320. Computers are well known in the art and willnot be described in detail herein.

The foregoing Detailed Description is to be understood as being in everyrespect illustrative and exemplary, but not restrictive, and the scopeof the invention disclosed herein is not to be determined from theDetailed Description, but rather from the claims as interpretedaccording to the full breadth permitted by the patent laws. It is to beunderstood that the embodiments shown and described herein are onlyillustrative of the principles of the present invention and that variousmodifications may be implemented by those skilled in the art withoutdeparting from the scope and spirit of the invention. Those skilled inthe art could implement various other feature combinations withoutdeparting from the scope and spirit of the invention.

1. A computer-implemented-method for rendering an image from a set ofk-space data, comprising the steps of: constructing a cost functionalbased at least in part on said set of k-space data, said cost functionalcomprisingλ∥φ(ƒ)∥₁ +∥W(Aƒ−y)∥₂ ² wherein ƒ is a vector representing said set ofimage-space data, λ is a Lagrange multiplier, y is a vector representingsaid set of k-space data, φ(ƒ) is an operator which transforms ƒ into asparse representation, A is a Fourier operator which maps said set ofimage-space data into said set of k-space data, and W is said weightingmatrix; iteratively solving for a set of image-space data, wherein saidstep of iteratively solving comprises the step of minimizing said costfunctional; and, rendering said image from said set of image-space data.2. The method of claim 1 wherein said step of iteratively solvingfurther comprises the steps of: generating a block-Toeplitz matrix basedat least in part on said weighting matrix; and, performing a pluralityof first iterations based at least in part on said block-Toeplitz matrixelements, wherein each first iteration comprises the steps of: (a)calculating the gradient of said cost functional; (b) performing aplurality of second iterations, wherein each second iteration comprisesthe steps of: (b1) performing a backtracking line search; and, (b2)calculating a value of said cost functional; and, (c) updating said setof image-space data.
 3. The method of claim 2, wherein saidblock-Toeplitz matrix comprises:A ^(†) DA wherein D=(AA^(†)) ⁻¹ .
 4. The method of claim 3 wherein saidgradient of said cost functional comprises:∇ψ=λ∇∥φ(ƒ)∥₁+2A ^(†) D(Aƒ−y).
 5. The method of claim 3 wherein said stepof iteratively solving further comprises the step of iteratively solvingthe equation: $\begin{matrix}{f_{k} \approx {f_{k - 1} + \left( {{{- A^{\dagger}}{DAf}_{k - 1}} + {A^{\dagger}{Dy}}} \right)}} \\{= {{A^{\dagger}{Dy}} + {\left( {I - {A^{\dagger}{DA}}} \right)f_{k - 1}}}}\end{matrix}$ wherein ƒ_(k) is the set of f-values after the k-thiteration and I is the identity matrix.
 6. The method of claim 1 wherein∥φ(ƒ)∥₁ is a total variation.
 7. An image processing system forrendering an image from a set of k-space data, comprising: means forconstructing a cost functional based at least in part on said set ofk-space data, said cost functional comprising:λ∥φ(ƒ)∥₁ +∥W(Aƒ−y)∥₂ ² wherein ƒ is a vector representing said set ofimage-space data, λ is a Lagrange multiplier, y is a vector representingsaid set of k-space data, φ(ƒ) is an operator which transforms f into asparse representation, A is a Fourier operator which maps said set ofimage-space data into said set of k-space data, and W is said weightingmatrix; means for iteratively solving for a set of image-space data,wherein said step of iteratively solving comprises the step ofminimizing said cost functional; and, means for rendering said imagefrom said set of image-space data.
 8. The image processing system ofclaim 7, further comprising: means for generating a block-Toeplitzmatrix based at least in part on said weighting matrix; and, means forperforming a plurality of first iterations based at least in part onsaid block-Toeplitz matrix elements, wherein each first iterationcomprises the steps of: (a) calculating the gradient of said costfunctional; (b) performing a plurality of second iterations, whereineach second iteration comprises the steps of: (b1) performing abacktracking line search; and, (b2) calculating a value of said costfunctional; and, (c) updating said set of image-space data.
 9. The imageprocessing system of claim 8 wherein said block-Toeplitz matrixcomprises:A ^(†) DA wherein D=(AA†) ⁻¹ .
 10. The image processing system of claim9 wherein said gradient of said cost functional comprises:∇ψ=λ∇∥φ(ƒ)∥₁+2A ^(†) D(Aƒ−y).
 11. The image processing system of claim9, further comprising means for solving the equation: $\;\begin{matrix}{f_{k} \approx {f_{k - 1} + \left( {{{- A^{\dagger}}{DAf}_{k - 1}} + {A^{\dagger}{Dy}}} \right)}} \\{= {{A^{\dagger}{Dy}} + {\left( {I - {A^{\dagger}{DA}}} \right)f_{k - 1}}}}\end{matrix}$ wherein ƒ_(k) is the set of f-values after the k-thiteration and I is the identity matrix.
 12. The image processing systemof claim 7 wherein ∥φ(ƒ)∥₁ is a total variation.
 13. A computer readablemedium storing computer program instructions for rendering an image froma set of k-space data, the computer program instructions defining thesteps of: constructing a cost functional based at least in part on saidset of k-space data, said cost functional comprisingλ∥φ(ƒ)∥₁ +∥W(Aƒ−y)∥₂ ² wherein ƒ is a vector representing said set ofimage-space data, λ is a Lagrange multiplier, y is a vector representingsaid set of k-space data, φ(ƒ) is an operator which transforms ƒ into asparse representation, A is a Fourier operator which maps said set ofimage-space data into said set of k-space data, and W is said weightingmatrix; iteratively solving for a set of image-space data, wherein saidstep of iteratively solving comprises the step of minimizing said costfunctional; and, rendering said image from said set of image-space data.14. The computer readable medium of claim 13 wherein said computerprogram instructions further comprise computer instructions defining thesteps of: generating a block-Toeplitz matrix based at least in part onsaid weighting matrix; and, performing a plurality of first iterationsbased at least in part on said block-Toeplitz matrix elements, whereineach first iteration comprises the steps of: (a) calculating thegradient of said cost functional; (b) performing a plurality of seconditerations, wherein each second iteration comprises the steps of: (b1)performing a backtracking line search; and, (b2) calculating a value ofsaid cost functional; (c) updating said set of image-space data.
 15. Thecomputer readable medium of claim 14, wherein said block-Toeplitz matrixcomprises:A ^(†) DA wherein D=(AA†) ⁻¹ .
 16. The computer readable medium of claim15, wherein said gradient of said cost functional comprises:∇ψ=λ∇∥φ(ƒ)∥₁+2A ^(†) D(Aƒ−y).
 17. The computer readable medium of claim15, wherein said computer instructions further comprise computerinstructions for solving the equationf_(k) ≈ f_(k − 1) + (−A^(†)DAf_(k − 1) + A^(†)Dy)   = A^(†)Dy + (I − A^(†)DA)f_(k − 1)wherein ƒ_(k) is the set of ƒ-values after the k-th iteration and I isthe identity matrix.
 18. The computer readable medium of claim 13wherein ∥φ(ƒ)∥₁ is a total variation.